#' Compute Krippendorff's alpha for multiple rating/coding/annotation data with
#'     nominal/categorical values/label classes.
#'     
#' Given a "reliability" data frame 
#'    (i.e., a matrix or data frame of coders-times-units recording coder j's rating of unit i),
#'    function computes Krippendorff's alpha.
#'    
#' @alpha a tabular data object 
#'     (2D-\code{\link{array}}, \code{\link{matrix}} or one inheriting from class \code{\link{data.frame}})
#'     representing coders in rows 1:J and items in columns 1:I.
#'     
#'     \bold{NOTE} that If a \code{data.frame} is provided, it is assumed that the first 
#'     column records coder ITs, however, so that items are represented in columns 2:(I+1).
#' 
#' @return A scalar numeric value (type \code{double}) equal to Krippendorf's alpha.
#' 
#' @note A more flexible implementation in Pyhton is available at \url{https://github.com/pln-fing-udelar/fast-krippendorff}
#' @source Example data from \url{https://en.wikipedia.org/wiki/Krippendorff%27s_alpha#A_computational_example}
#' 
#' @examples
#' \dontrun{
#'  reldat <- tibble::tribble(
#'    ~coder, ~u1, ~u2, ~u3, ~u4, ~u5, ~u6, ~u7, ~u8, ~u9, ~u10, ~u11, ~u12, ~u13, ~u14, ~u15,
#'       "A",  NA,  NA,  NA,  NA,  NA,  3L,  4L,  1L,  2L,   1L,   1L,   3L,   3L,   NA,   3L,
#'       "B",  1L,  NA,  2L,  1L,  3L,  3L,  4L,  3L,  NA,   NA,   NA,   NA,   NA,   NA,   NA,
#'       "C",  NA,  NA,  2L,  1L,  3L,  4L,  4L,  NA,  2L,   1L,   1L,   3L,   3L,   NA,   4L,
#'  )
#'  krippendorffs_alpha(reldat)
#' }
krippendorffs_alpha <- function(x) {
  
  if (!any(c(inherits(x, "array"), inherits(x, "matrix"), inherits(x, "data.frame"))))
    stop("`x` must be a 2D-array, matrix or inherit from 'data.frame'.")
  
  if (is.array(x) && length(dim(x)) != 2L)
    stop("`x` is an array but must be a 2-dimensional.")
    
  if (inherits(x, "data.frame"))
    x <- as.data.frame(x[-1])
  
  cats <- na.omit(unique(unlist(lapply(x, unique))))
  store <- matrix(0, nrow = length(cats), ncol = length(cats), dimnames = list(cats, cats))
  
  res <- lapply(
    x
    , function(
      .x # vector of coders' unit ratings
      , .store = store# initial coincidence matrix
    ) {
      # get non-missing ratings
      vals <- na.omit(.x)
      # get number of non-missing ratings
      m <- length(vals)
      # return if no pairable ratings
      if (m < 2L)
        return(.store)
      
      # get permutations of pairable ratings
      pairs <- gtools::permutations(length(vals), 2L, vals, set = FALSE)
      
      # iterate over rows and increment .store's elements in place
      apply(pairs, 1, function(.x) {.store[.x[1], .x[2]] <<- .store[.x[1], .x[2]] + 1L})
      
      # (element-wise) division by m-1 and return
      return(.store/(m-1))
  })
  
  # add units' concidence matrices element-wise 
  o <- Reduce("+", res)
  
  # compute observed disagreement 
  D_o <- sum(o[upper.tri(o)])
  
  # compute margins
  n_v <- rowSums(o)
  
  # compute expected disagreement 
  tmp <- combn(n_v, 2)
  D_e <- 1/(sum(n_v)-1)*sum(tmp[1, ]*tmp[2, ])
  
  # compute alpha
  alpha = 1-(D_o/D_e)
  
  # return
  return(alpha)
}

